ANL-HEP-CP-99-35 
NIKHEF/99-013 
ITP-SB-99-18 



Heavy Quark Production in Deep-Inelastic Scattering at 

HERA 

B.W. Harris , E. Laenen 6 , S. Moch 6 , J. Smith c 

a Argonne National Laboratory, 9700 S. Cass Ave., Argonne, IL 60439, USA 

b NIKHEF Theory Group, Kruislaan 409, 1098 SJ Amsterdam, The Netherlands 

c Institute for Theoretical Physics, S.U.N.Y. at Stony Brook, Stony Brook, NY 11794-3840, USA 

Abstract: We discuss two topics in the production of heavy quarks in deep- 
inelastic scattering: the next-to-leading order Monte-Carlo HVQDIS and the next- 
to-leading logarithmic resummation of soft gluon effects, including estimates of next- 
to-next-to-leading order corrections therefrom. 

1 Introduction 

Charm quarks produced in deep-inelastic scattering have been identified in sizable numbers by 
the HI and ZEUS collaborations at HERA [1], and considerably more charm (and bottom) 
data are anticipated. At the theoretical level the reaction has already been studied extensively. 
In the framework where the heavy quark is not treated as a parton, leading order (LO) [2] 
and next-to-leading order (NLO) [3] calculations of the inclusive structure functions exist. 
Moreover, LO (AROMA, RAPGAP) [4,5] and NLO (HVQDIS) [6,7] Monte-Carlo programs, 
allowing a much larger class of observables to be compared with data, have been constructed 
in recent years. The NLO QCD description agrees quite well with the HERA data. Most 
of the recent theoretical attention for this reaction has been in the context of variable flavor 
number schemes [8] , which we shall not address here. We shall review here two issues regarding 
heavy quark production in deep-inelastic scattering. In the next section we discuss two new 
features of HVQDIS, relevant to recent and future analyses; first, the inclusion of charmed- 
meson semi-leptonic decays, and second, a switch to describe (LO) D-D-jet final states. In the 
third section we review the methods and some key results of the next-to-leading logarithmic 
Sudakov resummation for DIS production of heavy quarks, and NNLO estimates derived from 
this resummation. 
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2 The NLO Monte-Carlo HVQDIS 



The program HVQDIS [6] is based on the next-to-leading order fully differential heavy quark 
contributions to the proton structure functions [7]. The basic components (in terms of virtual- 
photon-proton scattering) are the 2-to-2 body squared matrix elements through one-loop order 
and tree level 2-to-3 body squared matrix elements, for both photon-gluon and photon-light- 
quark initiated subprocesses. It is therefore possible to study single- and semi-inclusive pro- 
duction at next-to-leading order, and three body final states at leading order. The goal of the 
next-to-leading order calculation is to organize the soft and collinear singularity cancellations 
without loss of information in terms of observables that can be predicted. 

The subtraction method provides a mechanism for this cancellation. It allows one to isolate 
the soft and collinear singularities of the 2-to-3 body processes within the framework of dimen- 
sional regularization without calculating all the phase space integrals in a space-time dimension 
n^4. Expressions for the three-body squared matrix elements in the limit where an emitted 
gluon is soft appear in a factorized form where poles e = 2 — n/2 multiply leading order squared 
matrix elements. These soft singularities cancel upon addition of the interference of the leading 
order diagrams with the renormalized one-loop virtual diagrams. The remaining singularities 
are initial state collinear in origin wherein the three-body squared matrix elements appear in 
a factorized form, with poles in e multiply splitting functions convolved with leading order 
squared matrix elements. These collinear singularities are removed through mass factorization. 

The result of this calculation is an expression that is finite in four-dimensional space time. 
One can compute all phase space integrations using standard Monte-Carlo integration tech- 
niques. The final result is a program which returns parton kinematic configurations and their 
corresponding weights, accurate to 0(aa 2 s ). The user is free to histogram any set of infrared- 
safe observables and apply cuts, all in a single histogramming subroutine. Additionally, one 
may study heavy hadrons using the Peterson et al. model. Detailed physics results from this 
program are given in [9]. 

Below we discuss and give examples of two new options of HVQDIS version 1.3 1 : electrons 
from semileptonic decays of the charmed hadron, and a switch for retaining only three body 
final states. 

2.1 Semileptonic decays 

HVQDIS has been extended to include the electron from the semileptonic decay of the charmed 
hadron. In the rest frame of the decaying hadron, the electrons are distributed isotropically. 
Their momentum distribution is the weighted average of multiple decay modes of many different 
charmed hadrons, and has been deduced [10] from RAPGAP[5]. The parameterization is 
shown in fig. 1. The implementation used in HVQDIS1.3 (shown as points) was derived from 
the fit (dashed line) to the RAP GAP output (histogram). The overall normalization of the 
cross section is then fixed by the semileptonic branching ratio Br(c — > e) which we take to be 
9.5%. 

The inclusive transverse momentum and pseudo-rapidity distributions of the semileptonic 
decay electrons in the lab frame in deep inelastic scattering of 820 GeV protons with 27.5 GeV 
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Figure 1: Momentum distribution of semileptonic decay electrons from charmed hadrons pro- 
duced by RAPGAP in the rest frame of the hadron. Momenta distributed according to the fit 
are implemented in HVQDIS1.3. 

electrons in the kinematic range < y < 0.7 and 2 < Q 2 < 100 GeV 2 are shown in fig. 2. 
We also show the corresponding distributions for the parent parton (c-quark) and hadron 
(charmed- meson). The curves are produced using the next-to-leading order Gliick-Reya-Vogt 
1994 (GRV94) [11] parton distribution functions, a two-loop a s with n/ = 3 and Aq^^ = 248 
MeV, and m c = 1.35 GeV. The distributions of the charmed partons and hadrons are highly 
correlated because of the simple Peterson et al. fragmentation model. The semileptonic decay 
electrons are very soft, taking only a small portion of the hadron p t , and more central due to 
the isotropic nature of the decay. 



2.2 Three body final states 

For speed considerations it pays to add a switch to turn off all two body contributions (primarily 
the very slow virtual routines) when one is interested only in a manifestly three body observable. 
Such a switch has been added to HVQDIS1.3. We give here a sample of three body observables. 

The final state of interest is D-D-jet corresponding to the partonic states c-c-g and c-c-q. 
We begin by requiring the D, D and jet to be above some minimum transverse momentum 
(P t D > 1.2 GeV, Pp > 1.2 GeV, P t jet > 6 GeV), to be central (\r] D \ < 1.5, \r^\ < 1.5, 
\rjjet\ < 2.4), and to be isolated (R D jj > 0.7, Rojet > 0.7, Rjjj et > 0.7) in the lab frame. The 

cone size R^ = - Vj) 2 + (0i - 0j) 2 - 

The total cross section for the deep inelastic production of D-D-jet as a function of their 
invariant mass Mjjpj is shown in fig. 3 for < y < 0.7 and 2 < Q 2 < 100 GeV 2 . A one loop 

a s with rif = 3 and Aqq^ = 232 MeV, leading order GRV94[11] parton distributions, and 
m c = 1.35 GeV where used. The normalization of this LO curve has a large uncertainty as 
demonstrated by the various scale choices = {/U /2, /x , 2/i }, with /i = ^JQ 2 + m 2 + M^.. 
Also shown in the figure is a decomposition into the gluon and light-quark initiated subprocesses. 
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Figure 2: The inclusive transverse momentum (left) and psuedo-rapidity (right) distributions 
of the semileptonic electron, the parent parton (c-quark), and the hadron (charmed-meson) in 
the lab frame for < y < 0.7 and 2 < Q 2 < 100 GeV 2 . 




Figure 3: Total cross section for the deep inelastic production of D-D- jet as a function of 
their invariant mass M D -^- for < y < 0.7 and 2 < Q 2 < 100 GeV 2 . The central scale is 

Ho = ^JQ 2 + m 2 + M 2 D -^. . A decomposition into quark and gluon initiated subprocesses is also 
shown. 



The gluon initiated subprocess dominates due to the relatively large size of the gluon parton 
distribution function at small x. As another example, in the DDj center-of-mass frame we 
construct the Dalitz energy fractions Xi = 2Ei/M D -^-, (i = D, D, or j) that specify how much 
available energy is shared between the D,D, and jet. They satisfy xo + xjy + Xj = 2. The 
normalized cross section differential in io and Xj is shown in fig. 4. 
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Figure 4: Dependence of the DDj normalized cross section on the energy fraction x for the D 
hadron (left) and jet (right) for < y < 0.7 and 2 < Q 2 < 100 GeV 2 . 

3 Soft-gluon resummation 

As already remarked, existing NLO calculations for heavy quark electroproduction provide a 
solid theoretical perturbative QCD description [3,7,9] for this reaction, and agree well with 
present data [1]. At moderate Q 2 and x values larger than 0.01, the charm structure function 
-Fcharm is increasingly dominated by partonic processes near the charm quark pair production 
threshold. The large size of the gluon density f g (x,ii) for small momentum fractions x gives 
relatively large weight to such processes [12]. Although the QCD corrections at presently acces- 
sible x values are moderate (about 30-40%), with an increasing number of data to be gathered 
at higher x, it is worthwhile to have a closer look at such near-threshold subprocesses. In this 
kinematic region, the QCD corrections are dominated by large Sudakov double logarithms. Re- 
cently [13], these Sudakov logarithms have been resummed to all orders of perturbation theory, 
to next-to- leading logarithmic (NLL) accuracy, and, moreover, in single-particle inclusive (1PI) 
kinematics [14] 2 . Let us recall at this point the main results. First, the quality of the ap- 
proximation for the next- to-leading logarithmic threshold resummation was found to be clearly 
superior to leading logarithmic one. Furthermore, the resummation provided NNLO estimates 
[13], which were found to be sizable for x > 0.05. 

Below we give a synopsis of the analysis presented in [13]. We study electron proton scat- 
tering with the exchange of a single virtual photon, Q 2 = —q 2 , and a detected heavy quark (we 
concentrate on the charm quark case here) in the final state, i.e. the subprocess 

7 ( g ) + P(p) — Q(p 1 ) + X[Q], (1) 

where X denotes any additional hadrons, including the heavy anti-quark, in the final state 
and p\ = m 2 . The Mandelstam invariants, S' = (p + q) 2 + Q 2 , T"i = (p — pi) 2 — m 2 and 
Ui = (q — pi) 2 — m 2 can be used to define S4 = S' + 7\ + Ui, which vanishes at the hadronic 
threshold. The double differential heavy quark structure function dF 2 associated to (1) may be 

2 Analytical results for pair-inclusive kinematics are also given in [13] 
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written as 



d 2 F 2t p(x,S4,T 1 ,U 1 ,Q 2 ,m 2 ) _ I „ ) dz 2 /x s 4 t 1 u x Q 2 m 2 \ 

" ^ ™- 9 ^ ~ ' ^ W2,< ^' ^' ^' ^ ' 

(2) 

where z = —U\/{S' + Ti). The fi/p denote parton distributions in the proton at momentum 
fraction z and MS-mass factorization scale \i. The functions u 2ji describe the underlying hard 
parton scattering processes and depend on the partonic Mandelstam variables s',ti,ui and S4, 
which are derived from (1) by replacing the proton P by a parton of momentum k = zp. At n-th 
order in perturbation theory, the gluonic hard part u 2t9 in eq.(2) typically depends on singular 
distributions a"[ln 2n ~ 1 ~ fc (s 4 /m 2 )/s4] + , k = 0,1,..., that must be resummed. Contributions 
from light initial-quark states are neglected, as they are about 5% at NLO. 

The resummation of threshold logarithms rests upon the factorization of the kinematics 
and dynamics of the relevant degrees of freedom near threshold [15,16]. The dynamical factors 
involved can be each be assigned a kinematic weight Wi that is defined to vanish at threshold. 
For dF 2) p in eq.(2), the factorization of the kinematics implies that these weights sum to the 
overall inelasticity near threshold: S^/m 2 ~ wi(— u\)/m 2 + w s , with wi — 1 — z and w s = 
s^/m 2 . Correspondingly, the infrared regulated partonic structure function dF 2t9 factorizes 
into functions that individually organize contributions from these near-threshold degrees of 
freedom. Thus, there is here a function i]) g / g that sums the singular distributions from incoming 
collinear gluons, and a soft function S that organizes those due to soft gluons not collinear to 
the incoming parton. Finally, there is a hard function H 2g incorporating only regular short- 
distances corrections. Replacing the proton in eq. (2) by a gluon, and passing to Laplace- moment 
space, f(N) = f£°dwexp[—Nw}f(w), this gives [13] 



I N — — 0— — \ H ( — — ^— — 



2 m 2 N 



ti u x Q 2 



lpg/g(N u ,p- C/fi) 
4>g/g{N u ,fi) 



m 



s {w ( ' (3) 



where (fi g / g is the usual MS-distribution from mass factorization and N u = N(—u\)/m 2 . In 
moment space, the Sudakov logarithms appear as factors a n s ln 2n "W, with i = 0, 1 for NLL 
accuracy. The iV-dependence in eq.(3) exponentiates for each function individually. All leading 
logarithms (LL) are exclusively contained in ip g / g , which is a gluon distribution at fixed fraction 
of p ■ ( and can be defined as an operator matrix element. It depends on a time-like vector 
( = p 2 /m {p 2 is the heavy antiquark momentum). Its collinear poles are canceled by 4>g/g- The 
threshold logarithms in ip g / g are resummed in analogy to the Drell-Yan process [15], while all 
scale dependence of ip g / g and (f> g / g is governed by renormalization group equations (RGE) with 
anomalous dimensions 7^ = 2 r y g and 7 s / 9 [17,13]. 

The soft function S requires renormalization, since it is defined as a composite operator, 
that connects Wilson lines in the direction of the scattering partons [17,18,19]. Its RGE, 
n{d/dn)\nS(N) = — 2Rer 5 , resums all threshold logarithms in S. Its gauge dependence 
cancels precisely that of ip g / g - The soft anomalous dimension T s is to order a s 



r s (e* a ) = 




In 



Am 



(4) 



with = ^l-Am 2 /s and L p = (1 - 2 m 2 /s){ln(l - /3)/(l + (5) + m}/(3. 
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The final result for the hard scattering function oj 2 , g in moment space resums all large 
logarithms in single-particle inclusive kinematics up to NLL accuracy. Combining the resummed 
■ipg/g with the integrated RGE for S, we obtain for u)2, g [13] 



~ In — — 91 — 

W2 ' g { V V /i 2 ' /i 2 



(5) 



H. 



ti ui Q 2 



TYl 

\ m z m z m 2 y jy 



exp < — 2 



%-7 fl (a a (A*')) 
A* 



x exp 



d A 1 

yA (9) (a s (v / Am)) + -z/ (g) (a s (u?m)) 



x 



{/"(— 

^ 

exp | J -^2Rer s (a s (pi')) Lxp j -2 j ( 7fl (a s (//)) - lg/g( N u, « s (a*'))) 

>■ m . ^ s ^ it. ^ 



The second exponent gives the leading iV-dependence of the ratio ipg/g/4>g/ g with z/( 9 )(o; s ) = 
2C A a 8 /n, A to) (a a ) = C A (a s /vr) + (C A ^/2)(a s /7r) 2 and K = C A (67/18 - tt 2 /6) - 5/9n f [20]. 
For NLL Sudakov resummation, the product H 2g • S on the first line of eq. (5) is determined 
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Figure 5: (a): The ^-dependence of the coefficient functions c 2g (77, £)> k = 0, 1, 2 for Q 2 = 10 GeV 2 

(k 0) 

with m = 1.5 GeV. Plotted are the exact results for c 2g , k = 0, 1 (solid lines), the LL approximation 

to c^'g ^ (dotted line) the NLL approximation to (dashed line), the LL approximation to c^ "* 

(dash-dotted line) and the NLL approximation to (long dashed line), (b): F2 harm (x, Q 2 , z max ) as 
a function of z max at NLO with the CTEQ4M gluon PDF, x = 0.01, m = 1.6 GeV, Q 2 = 10 GeV and 
H = m (upper three curves), fx = VQ 2 + 4m 2 (1 ower three curves), rescaled by a factor of 0.5. Plotted 
are: The exact results (solid lines), the LL approximations (dotted lines) and the NLL approximations 
(dashed lines). 



The resummed result for u) 2iff in eq. (5) may be used as a generating functional for fixed 
order approximate perturbation theory by re-expanding u 2 ^ g to NLO and NNLO and inverting 
the Laplace transform. After insertion of eq. (5) into eq. (2) and integration over Ti,Ui, we 
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Figure 6: (a): The x-dependence of the ratios F&^JFgfffy (solid line) and F% h ™™ LO) /F* h ™™ 
(dashed line) with F 2 ch ( ^™ LO) in the improved NLL approximation (exact NLO result plus NLL ap- 
proximate NNLO result) with CTEQ4M gluon PDF, /j = m = 1.6 GeV and Q 2 = lOGeV. (b): The 
differential distribution dF£ harm / dp T as a function of p T with the CTEQ4M gluon PDF, x = 0.01, 
m = 1.6 GeV, Q 2 = 10 GeV and scale choice fi = ^Q 2 + 4(m 2 + p 2 -). Plotted are: The exact re- 
sult (solid line) at NLO, the LL approximation at NLO (dotted line), the NLL approximation at 
NLO (dashed line) and at NNLO the improved NLL approximation (exact NLO result plus NLL 
approximate NNLO result) (long dash dotted line). 



may expand the structure function as 

Ft* m (x,Q 2 ,m 2 ) = a f\ e "f j dzf g/P (z^ 2 )±(Ana s (pt)) k j:c^\r l ,0ln 1 ^ (6) 

47r m i x k=o i=o m 

where a = (Q 2 + Am 2 ) /Q 2 and e c = 2/3. 

The quality of the NLL approximation eq.(5) can then be investigated by comparing exact, 
LL and NLL approximation to the gluon coefficient functions c^J\ which is done in fig. 5(a). 
The functions depend on the scaling variables 

S 1 * Q2 (7\ 

Am z m z 

where r\ is a direct measure of the distance to the partonic threshold. 

Fig. 5(a) reveals that, although at one loop the LL accuracy provides a good approximation 
for very small 77, the NLL approximation is excellent over a much wider range in 77, up to values 
of about 10 (the same actually holds true for the c£J\ I > 0, k < 2). We also show c^ "*, which 

has more structure than in the curves. 

To address the threshold sensitivity of the integrated charm structure function to threshold 
processes, we perform the integral over z in eq. (2) only up to a value z max , and plot the integral 
then as a function of z max . In this way we can see where the integral eq. (2) acquires most of its 
value. The physical structure function corresponds to z max = 1. Fig. 5(b) demonstrates that 
for x = 0.01 _P 2 charm is mostly determined by partonic processes close to threshold. 
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In fig. 6a we display at a fixed value of the factorization scale /j, — m over a range of x, 
0.003 < x < 0.3 the effect of the NNLO corrections. We plot the K-factors F 2 ch ( ^ i0) /^(tvlo) 
and, for comparison, also F^n^o)/ ^^(lo) 3 - At NNLO we have taken the improved NLL 
approximation to ,F 2 charm (the exact NLO result plus the NLL approximate NNLO result). We 
see that particularly for smaller x, the size of the NNLO corrections is negligible, the K-factor 
being close to one, whereas for larger x, their overall size is still quite big, almost a factor of 2 
at x — 0.1. 

Finally, in fig. 6b we show the NLO results as a function of pr for x = 0.01, m = 1.6 GeV, 
Q 2 = 10 GeV and fj, = \JQ 2 + 4(m 2 + pj). At NLO, we compare our LL and NLL approximate 
results with the exact results of the second Ref. in [3]. We see that the exact curves are 
reproduced well both in shape and magnitude by our NLL approximations, whereas the curves 
for LL accuracy systematically underestimate the true result. We also display the improved 
NLL approximation to <iF 2 charm / 'dpx at NNLO, which contains sizeable contributions to the 
value of the maximum increases by 40% - 50%. 



4 Conclusions 



Driven by the ever increasing variety and quantity of deep-inelastic charm production data 
from the HI and ZEUS experiments at HERA, we have updated and reviewed two important 
tools: soft gluon threshold resummation and the next-to-leading order Monte-Carlo HVQDIS. 
The addition of semileptonic decay electron information and an option for only three body final 
states to the Monte-Carlo will enhance future physics analysis options. Soft gluon threshold 
resummation, on the other hand, teaches us about the size of the terms neglected in fixed order 
calculations. 
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